################################################################################
#####15. Multilevel Models and Marginal Effect Plots############################
################################################################################

#####Load data from MA2_Data_Generation#####

pols <- readRDS('pols_bjps')

#Load nlme package#

library(nlme)

#####Main Models#####

#####No Interactions#####

no.int1e <- lme(scale_econsdw ~ scale_ggini_inc1 +
                         scale_socioeconomicsal +
                         scale_dispgini +
                         scale_enep +
                         scale_galtansdw +
                  scale_mag +
                  scale_unemp,
                       random = ~ 1 |  eastwest/ country.name,
                       data = pols,
                       na.action = na.omit,
                method = "ML")

no.int2e <- lme(scale_econsdw ~ scale_ggini_inc2 +
                 scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
                  scale_mag +
                  scale_unemp,
               random = ~ 1 | eastwest / country.name , 
               data = pols,
               na.action = na.omit,
               method = 'ML')

no.int3e <- lme(scale_econsdw ~ scale_ggini_inc3 +
                 scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
                  scale_mag +
                  scale_unemp,
               random = ~ 1 | eastwest / country.name , 
               data = pols,
               na.action = na.omit,
               method = 'ML')

no.int4e <- lme(scale_econsdw ~ scale_ggini_inc4 +
                 scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
                  scale_mag +
                  scale_unemp,
               random = ~ 1 |  eastwest / country.name ,
               data = pols,
               na.action = na.omit,
               method = 'ML')

no.int5e <- lme(scale_econsdw ~ scale_gcov_inc5 +
                 scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
                  scale_mag +
                  scale_unemp,
               random = ~ 1 |  eastwest / country.name , 
               data = pols,
               na.action = na.omit,
               method = 'ML')

#####All Interactions#####

all1e <- lme(scale_econsdw ~ scale_ggini_inc1 +
                 scale_socioeconomicsal +
                 scale_ggini_inc1 *scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
               scale_mag+
               scale_unemp +
                 scale_ggini_inc1 *scale_dispgini +
                 scale_socioeconomicsal * scale_dispgini +
                 scale_ggini_inc1 * scale_dispgini * scale_socioeconomicsal,
               random = ~ 1 | eastwest / country.name ,
               data = pols,
               na.action = na.omit,
             
             method = 'ML')

all2e <- lme(scale_econsdw ~ scale_ggini_inc2 +
                 scale_socioeconomicsal +
                 scale_ggini_inc2 * scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
               scale_mag+
               scale_unemp +
                 scale_ggini_inc2 * scale_dispgini +
                 scale_socioeconomicsal * scale_dispgini +
                 scale_ggini_inc2 * scale_dispgini * scale_socioeconomicsal,
               random = ~ 1 | eastwest / country.name , 
               data = pols,
               na.action = na.omit,
             method = 'ML')

all3e <- lme(scale_econsdw ~ scale_ggini_inc3 +
                 scale_socioeconomicsal +
                 scale_ggini_inc3 * scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
               scale_mag+
               scale_unemp +
                 scale_ggini_inc3 * scale_dispgini +
                 scale_socioeconomicsal * scale_dispgini +
                 scale_ggini_inc3 * scale_dispgini * scale_socioeconomicsal,
               random = ~ 1 | eastwest / country.name , 
               data = pols,
               na.action = na.omit,
             method = 'ML')

all4e <- lme(scale_econsdw ~ scale_ggini_inc4 +
                 scale_socioeconomicsal +
                 scale_ggini_inc4 * scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
               scale_mag+
               scale_unemp+
                 scale_ggini_inc4 * scale_dispgini +
                 scale_socioeconomicsal * scale_dispgini +
                 scale_ggini_inc4 * scale_dispgini * scale_socioeconomicsal,
               random = ~ 1 | eastwest / country.name , 
               data = pols,
               na.action = na.omit,
             method = 'ML')

all5e <- lme(scale_econsdw ~ scale_ggini_inc5 +
                 scale_socioeconomicsal +
                 scale_ggini_inc5 * scale_socioeconomicsal +
                 scale_dispgini +
                 scale_enep +
                 scale_galtansdw +
               scale_mag+
               scale_unemp +
                 scale_ggini_inc5 * scale_dispgini +
                 scale_socioeconomicsal * scale_dispgini +
                 scale_ggini_inc5 * scale_dispgini * scale_socioeconomicsal,
               random = ~ 1 | eastwest / country.name , 
               data = pols,
               na.action = na.omit,
             method = 'ML')


#####Table Output#####
#install.packages('stargazer')
library(stargazer)
#install.packages('MuMIn')
library(MuMIn)
#install.packages('dplyr')
library(dplyr)
#install.packages('ggplot2')
library(ggplot2)
#install.packages('ggrepel')
library(ggrepel)
#install.packages('MASS')
library(MASS)
#install.packages('gridExtra')
library(gridExtra)

#extract R2 for all plausible sets of values#

r1o <- r.squaredGLMM(no.int1e)
r2o <- r.squaredGLMM(no.int2e)
r3o <- r.squaredGLMM(no.int3e)
r4o <- r.squaredGLMM(no.int4e)
r5o <- r.squaredGLMM(no.int5e)

#bind them together#

r.o <- rbind(r1o, r2o, r3o, r4o, r5o)

#Take mean to get marginal and conditional R2$

mean(r.o[,1])
mean(r.o[,2])

#Essentially repeat for model fit statistics#

fit1 <- c(summary(no.int1e)$logLik, summary(no.int1e)$AIC, summary(no.int1e)$BIC)
fit2 <- c(summary(no.int2e)$logLik, summary(no.int2e)$AIC, summary(no.int2e)$BIC)
fit3 <- c(summary(no.int3e)$logLik, summary(no.int3e)$AIC, summary(no.int3e)$BIC)
fit4 <- c(summary(no.int4e)$logLik, summary(no.int4e)$AIC, summary(no.int4e)$BIC)
fit5 <- c(summary(no.int5e)$logLik, summary(no.int5e)$AIC, summary(no.int5e)$BIC)

fit.stats <- rbind(fit1, fit2, fit3, fit4, fit5)

print(mean(fit.stats[,1])) #log liklihood
print(mean(fit.stats[,2])) #AIC
print(mean(fit.stats[,3])) #BIC

#Standard Deviation of random intercepts - taken from summary for each model#
print(mean(c(0.7509236,0.7483622, 0.7493955, 0.750254,  0.7447141 ))) #Country
print(mean(c(0.0001075551,9.489712e-05, 8.388745e-05,8.509325e-05, 0.0001142614))) #Region



##Regression Table and Dot Plot##

#Extract predcitor info from each model#

m1 <- as.data.frame(summary(no.int1e)$tTable) 

m2 <- as.data.frame(summary(no.int2e)$tTable) %>%
  dplyr::rename(Value2 = Value, 
                Std.Error2 = Std.Error, 
                DF2 = DF, 
                `t-value2` = `t-value`, 
                `p-value2` = `p-value`)
m3 <- as.data.frame(summary(no.int3e)$tTable) %>%
  dplyr::rename(Value3 = Value, 
                Std.Error3 = Std.Error, 
                DF3 = DF, 
                `t-value3` = `t-value`, 
                `p-value3` = `p-value`)
m4 <- as.data.frame(summary(no.int4e)$tTable) %>%
  dplyr::rename(Value4 = Value, 
                Std.Error4 = Std.Error, 
                DF4 = DF, 
                `t-value4` = `t-value`, 
                `p-value4` = `p-value`)
m5 <- as.data.frame(summary(no.int5e)$tTable) %>%
  dplyr::rename(Value5 = Value, 
                Std.Error5 = Std.Error, 
                DF5 = DF, 
                `t-value5` = `t-value`, 
                `p-value5` = `p-value`)

#Bind data together#

dot.data <- cbind(m1,m2, m3,m4,m5)

#Reorder the data from to be easier ot read#

dot.data <- dot.data[ , order(names(dot.data))]

#Combine data to create pooled model#

dot.table <- dot.data %>%
  transmute(pe = rowMeans(dplyr::select(dot.data, Value, Value2, Value3, Value4, Value5)),
            se = rowMeans(dplyr::select(dot.data, Std.Error, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            DF = rowMeans(dplyr::select(dot.data, DF, DF2, DF3, DF4, DF5)),
            t.value = rowMeans(dplyr::select(dot.data, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data, `p-value`, `p-value2`, `p-value3`, `p-value4`, `p-value5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

#Create Predictor Column

dot.table$Predictor <- c("Intercept", "Income.Dif", "Socio-Economic Salience", "Disposable Income Inequality (Gini)", "ENEP", "GALTAN Polarization", "District Magnitude", "Unemployment")

#Create separate version of making dotplot#

dot.plot <- dot.table

#Reorder predictors#

dot.plot$Predictor <- factor(dot.plot$Predictor, 
                             levels = c("Income.Dif", "Socio-Economic Salience", "Disposable Income Inequality (Gini)", "GALTAN Polarization", "ENEP", "District Magnitude", "Unemployment", "Intercept"))

#reverse the direction of the predictors#

dot.plot$Predictor <- factor(dot.plot$Predictor, levels=rev(levels(dot.plot$Predictor)))

dot.plot.lm <- dot.plot[c(2:4),]

#Make the dotplot for key variables in the no interaction model - Not in the Paper #

dot.plot.lm %>% ggplot(aes(x = Predictor, y = pe)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
  geom_linerange(aes(ymin = lwr, ymax = up))+
  coord_flip() +
  labs(title = "Point Estimates of Model Regressing Predictors on Economic Polarization (n = 82)", y = "Point Estimates (95% Confidence Interval)", x = "Predictor")+
  theme_bw()

rownames(dot.table) <- dot.table$Predictor
dot.table <- dot.table[,-8]

#Regression Results for Model 3 in Table 1#
print(xtable::xtable(dot.table, digits = 3))

#####Point Estimates with all Interactions#####
r1o <- r.squaredGLMM(all1e)
r2o <- r.squaredGLMM(all2e)
r3o <- r.squaredGLMM(all3e)
r4o <- r.squaredGLMM(all4e)
r5o <- r.squaredGLMM(all5e)

r.o <- rbind(r1o, r2o, r3o, r4o, r5o)
mean(r.o[,1])
mean(r.o[,2])

#Standard deviations of random intercepts - from model summaries#

print(mean(c(5.48395e-05,6.989997e-05,6.111942e-05,5.847572e-05,7.037325e-05 )))#Region
print(mean(c(0.8517018,0.8368459,0.8440248,0.8469762,  0.8111853 )))#country


#Essentially repeat for model fit statistics#

fit1a <- c(summary(all1e)$logLik, summary(all1e)$AIC, summary(all1e)$BIC)
fit2a <- c(summary(all2e)$logLik, summary(all2e)$AIC, summary(all2e)$BIC)
fit3a <- c(summary(all3e)$logLik, summary(all3e)$AIC, summary(all3e)$BIC)
fit4a <- c(summary(all4e)$logLik, summary(all4e)$AIC, summary(all4e)$BIC)
fit5a <- c(summary(all5e)$logLik, summary(all5e)$AIC, summary(all5e)$BIC)

fit.statsa <- rbind(fit1a, fit2a, fit3a, fit4a, fit5a)

print(mean(fit.statsa[,1])) #log liklihood
print(mean(fit.statsa[,2])) #AIC
print(mean(fit.statsa[,3])) #BIC


##Regression Table and Dot Plot##

m1o <- as.data.frame(summary(all1e)$tTable) 

m2o <- as.data.frame(summary(all2e)$tTable) %>%
  dplyr::rename(Value2 = Value, 
                Std.Error2 = Std.Error, 
                DF2 = DF, 
                `t-value2` = `t-value`, 
                `p-value2` = `p-value`)
m3o <- as.data.frame(summary(all3e)$tTable) %>%
  dplyr::rename(Value3 = Value, 
                Std.Error3 = Std.Error, 
                DF3 = DF, 
                `t-value3` = `t-value`, 
                `p-value3` = `p-value`)
m4o <- as.data.frame(summary(all4e)$tTable) %>%
  dplyr::rename(Value4 = Value, 
                Std.Error4 = Std.Error, 
                DF4 = DF, 
                `t-value4` = `t-value`, 
                `p-value4` = `p-value`)
m5o <- as.data.frame(summary(all5e)$tTable) %>%
  dplyr::rename(Value5 = Value, 
                Std.Error5 = Std.Error, 
                DF5 = DF, 
                `t-value5` = `t-value`, 
                `p-value5` = `p-value`)

dot.data.o <- cbind(m1o,m2o, m3o,m4o,m5o)

dot.data.o <- dot.data.o[ , order(names(dot.data.o))]

dot.table.o <- dot.data.o %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.o, Value, Value2, Value3, Value4, Value5)),
            se = rowMeans(dplyr::select(dot.data.o, Std.Error, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            DF = rowMeans(dplyr::select(dot.data.o, DF, DF2, DF3, DF4, DF5)),
            t.value = rowMeans(dplyr::select(dot.data.o, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.o, `p-value`, `p-value2`, `p-value3`, `p-value4`, `p-value5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.o$Predictor <- c("Intercept", "Income.Dif", "Socio-Economic Salience", "Disposable Income Inequality (Gini)", "ENEP", "GALTAN Polarization", "District Magnitude", "Unemployment", "Income.Dif:Salience", "Income.Dif:Gini", "Salience:Gini", "Income.Dif:Salience:Gini")

dot.plot.o <- dot.table.o

dot.plot.o$Predictor <- factor(dot.plot.o$Predictor, 
                               levels = c( "Income.Dif", "Socio-Economic Salience", "Disposable Income Inequality (Gini)",  "GALTAN Polarization","ENEP", "District Magnitude", "Unemployment", "Income.Dif:Salience", "Income.Dif:Gini", "Salience:Gini", "Income.Dif:Salience:Gini", "Intercept"))

dot.plot.o$Predictor <- factor(dot.plot.o$Predictor, levels=rev(levels(dot.plot.o$Predictor)))
dot.plot.lm2 <- dot.plot.o[c(2:4, 9:12),]

coef_plot <- rbind(dot.plot.lm, dot.plot.lm2)
coef_plot$model <- c('No Interactions', 'No Interactions', 'No Interactions', rep('Interactions', 7))

coef_plot$Predictor <- factor(coef_plot$Predictor, 
                               levels = rev(c( "Disposable Income Inequality (Gini)", "Socio-Economic Salience", "Income.Dif", "Income.Dif:Gini", "Salience:Gini","Income.Dif:Salience", "Income.Dif:Salience:Gini")))

coef_plot$model <- factor(coef_plot$model, 
                          levels = c("No Interactions", "Interactions"))

coef_plot%>% ggplot(aes(x = Predictor, y = pe)) + 
  geom_point(color = 'black') + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed', color = 'black') +
  geom_linerange(aes(ymin = lwr, ymax = up), color = 'black')+
  coord_flip() +
  labs(title = "Point Estimates from Multi-Level Models 3 and 4 (Table 1 in Paper)", y = "Point Estimates (95% Confidence Interval)", x = "")+
  theme_bw()+
  facet_wrap( ~ model)

rownames(dot.table.o) <- dot.table.o$Predictor
dot.table.o <- dot.table.o[,-8]

#Regression results for Model 4 in Table 1#
print(xtable::xtable(dot.table.o, digits = 3))

#####MLM Interaction Plots#####
#extract fixed effects and variance-covariance matrix from model#
mu <- fixef(all1e)
sigma <- vcov(all1e)

#simulate 500 draws from a multivariate random normal distribution based on the model#
test <- as.data.frame(mvrnorm(500, mu = mu, Sigma = sigma))

######Figure 2 - Marginal Effect of Disposable Income Inequality on Party Polarization######

pols_rug <- readRDS('pols_bjps') %>%
  dplyr::select(scale_econsdw,scale_dispgini, 
                scale_ggini_inc1,
                scale_socioeconomicsal,
                scale_galtansdw ,
                scale_enep ,
                scale_unemp) %>%
  na.omit

nums <- seq(from = min(pols$scale_ggini_inc1, na.rm = T), to = max(pols$scale_ggini_inc1, na.rm = T), length.out = 500)

marg <- data.frame(x = nums)

for(i in 1:500){
  for(j in 1:500){
    marg[i,j+1] <- test[j,4] + test[j,10]*marg[i,1] + test[j,11]*-1 + test[j,12]*marg[i,1]*-1
  }
}

econsalw <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:500){
  econsalw[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:501], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:501], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:501], probs = .975))
}

marg.plot.dat <- data.frame(econsalw = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

lwr.plot.v <- ggplot(marg.plot.dat, aes(x = econsalw, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = -1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

#Middle Plot#

for(i in 1:500){
  for(j in 1:500){
    marg[i,j+1] <- test[j,4] + test[j,10]*marg[i,1] + test[j,11]*0 + test[j,12]*marg[i,1]*0
  }
}

econsalw <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:500){
  econsalw[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:501], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:501], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:501], probs = .975))
}

marg.plot.dat <- data.frame(econsalw = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

mid.plot.v <- ggplot(marg.plot.dat, aes(x = econsalw, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = 0", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

#High plot#

for(i in 1:500){
  for(j in 1:500){
    marg[i,j+1] <- test[j,4] + test[j,10]*marg[i,1] + test[j,11]*1 + test[j,12]*marg[i,1]*1
  }
}

econsalw <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:500){
  econsalw[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:501], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:501], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:501], probs = .975))
}

marg.plot.dat <- data.frame(econsalw = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

upr.plot.v <- ggplot(marg.plot.dat, aes(x = econsalw, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = 1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

#Combine Plots#

mlm_full_sample_sortingx<- grid.arrange(lwr.plot.v, mid.plot.v, upr.plot.v, nrow = 1, top = "", left = "Marginal Effect of Disposable Income Inequality", bottom = "Sorting by Income")

#ggsave("figure_A3.pdf", plot = mlm_full_sample_sortingx, width = 10, height = 4, units = "in")

#####Figure 3 - Marginal Effect of Disposable Income Inequality by Economic Salience######

nums <- seq(from = min(pols$scale_socioeconomicsal, na.rm = T), to = max(pols$scale_socioeconomicsal, na.rm = T), length.out = 500)

marg <- data.frame(x = nums)

for(i in 1:500){
  for(j in 1:500){
    marg[i,j+1] <- test[j,4] + test[j,11]*marg[i,1] + test[j,10]*-1 + test[j,12]*marg[i,1]*-1
  }
}

econsalw <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:500){
  econsalw[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:501], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:501], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:501], probs = .975))
}

marg.plot.dat <- data.frame(econsalw = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

lwr.plot.v2 <- ggplot(marg.plot.dat, aes(x = econsalw, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Party Income Differentiation = -1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

#Middle Plot#

for(i in 1:500){
  for(j in 1:500){
    marg[i,j+1] <- test[j,4] + test[j,11]*marg[i,1] + test[j,10]*0 + test[j,12]*marg[i,1]*0
  }
}

econsalw <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:500){
  econsalw[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:501], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:501], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:501], probs = .975))
}

marg.plot.dat <- data.frame(econsalw = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

mid.plot.v2 <- ggplot(marg.plot.dat, aes(x = econsalw, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Party Income Differentiation = 0", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

#High plot#

for(i in 1:500){
  for(j in 1:500){
    marg[i,j+1] <- test[j,4] + test[j,11]*marg[i,1] + test[j,10]*1 + test[j,12]*marg[i,1]*1
  }
}

econsalw <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:500){
  econsalw[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:501], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:501], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:501], probs = .975))
}

marg.plot.dat <- data.frame(econsalw = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

upr.plot.v2 <- ggplot(marg.plot.dat, aes(x = econsalw, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Party Income Differentiation = 1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

#Combine Plots#

mlm_full_sample_salx <- grid.arrange(lwr.plot.v2, mid.plot.v2, upr.plot.v2, nrow = 1, top = "", left = "Marginal Effect of Disposable Income Inequality", bottom = "Salience of Economic Issues")

#ggsave("figure_A4.pdf", plot = mlm_full_sample_salx, width = 10, height = 4, unit = "in")
